Evolution on a Rugged Landscape: Pinning and Aging 
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Population dynamics on a rugged landscape is studied analytically and numerically within a simple 
discrete model for evolution of A'' individuals in one-dimensional fitness space. We reduce the set of 
master equations to a single Fokker-Plank equation which allows us to describe the dynamics of the 
population in terms of thermo-activated Langevin diffusion of a single particle in a specific random 
potential. We found that the randomness in the mutation rate leads to pinning of the population 
and on average to a logarithmic slowdown of the evolution, resembling aging phenomenon in spin 
glass systems. In contrast, the randomness in the replication rate turns out to be irrelevant for 
evolution in the long-time limit as it is smoothed out by increasing "evolution temperature". The 
analytic results are in a good agreement with numerical simulations. 

PACS: 87.10.+e,82.20.Mj,64.60.Cn,05.40.+,i 



Recently, a simple model was introduced to describe 
the evolution of a finite population of mutating species 
in a one-dimensional fitness space \V^. Every species is 
characterized by a fitness variable which controls its repli- 
cation rate, and mutations move the species equiproba- 
bly up or down along this fitness axis. Analysis of this 
model leads to a two-staged dynamics of a population 
initially spread within some fitness range. At the first 
(fast) stage, the population forms a universal pulse-like 
distribution in the fitness space, which usually is accom- 
panied by the fast growth of the mean fitness. At the 
second stage, the pulse propagates toward higher fitness 
due to mutations, new more fit mutants are constantly 
generated, and less fit die due to the constant popula- 
tion size constraint. The mean fitness grows linearly in 
time. The non-trivial scaling of the fitness growth with 
the mutation strength and the population size was found 
in Refs. |l|J^ and appears to be in a qualitative agreement 
with the evolution dynamics of RNA virus [^ . 

Clearly, this model contains an implicit assumption 
that underlying fitness landscape in the configurational 
space is smooth. In reality, fitness landscapes likely ex- 
hibit substantial degree of ruggedness [Q, i.e. a proba- 
bility to find a more (or less) fit mutant varies for dif- 
ferent genomes. Moreover, for rugged landscapes there 
may be genomes which are more fit than any of their 
nearest neighbors in the sequence space which can be 
accessed as a result of a point mutation. Thus, there 
are local maxima of fitness. Rugged fitness landscapes 
are usually considered in multi-dimensional configura- 
tional spaces of individual species where the fitness is 
a complicated function of the detailed structure of the 
genome. A popular A'^if -model for the relation between 
genomic structure and fitness was introduced by Kauff- 
man H . Many other more complicated or realistic mod- 
els of rugged landscapes are described in the literature 
(see, e.g., ffl). 



Our model of one-dimensional fitness space is essen- 
tially different from these models since individual species 
with different internal structure but identical fitness are 
considered indistinguishable and therefore assigned to 
the same location along the fitness axis [Q. This model 
obviously cannot account for the evolution of the internal 
structure in the population, in particular, the famous er- 
ror catastrophe which occurs at high mutation rate and 
leads to accumulation of "bad mutants" and perpetual 
loss of heritable genetic information m . Nevertheless our 
model of evolution in a one-dimensional fitness space is 
capable of reproducing important dynamic features of 
evolution found in more realistic models. Due to its rel- 
ative simplicity, it is more amenable to analytical treat- 
ment then multi-dimensional configurational space mod- 
els. 

The goal of this paper is to study the dynamics of the 
population evolution in the framework of the model [y||] 
which however takes into account ruggedness of the un- 
derlying fitness landscape. Random fiuctuations which 
are considered to be functions of one fitness variable char- 
acterizing the genotype as a whole, modify locally both 
replication rate as well as mutation probability. We as- 
sume that the replication rate is a sum of linear func- 
tion of fitness and random fiuctuations. Linear part pro- 
vides the selective pressure driving the population to- 
wards higher fitness, as in smooth landscape case. Ran- 
dom fluctuations however provide trapping of the popula- 
tion near local maxima of replication rate and on average 
slow down the mean fitness growth. Probability of mu- 
tation up and down in fitness at a given state also has 
fluctuating part. We assume that these fluctuations are 
statistically independent of fluctuations of the replication 
rate. Based on master equations for underlying Markov 
process, we derive the Fokker-Plank equation for proba- 
bility distribution of species in the population. We find 
that initial growth of fitness is dominated by the fiuc- 



tuations of the replication rate, but at long-time limit 
these fluctuations become irrelevant. The long-time limit 
is completely determined by the fluctuations of the mu- 
tation rate which lead to logarithmic slowdown of the 
fitness growth. This phenomenon is analogous to the 
thermally-activated creep of a particle in quenched ran- 
dom (pinning) potential [^. 

Let us first specify the model (cf. Rcf. ||^). Consider a 
population of N individuals which can replicate accord- 
ing to their replication rate R, and babies can mutate, 
thus changing their fitness relative to parents. We assume 
asexual replication, i.e. any individual has a chance to re- 
produce independently. Once a baby is born, some mem- 
ber of the whole population (including the new baby) 
which is picked at random, is eliminated to preserve the 
constant population size. We assume that the replication 
rate i? is a function of only one independent variable x 
(fitness), which has a linearly growing component and 
a random component, R{x) — x + ^r{x) (without loss of 
generality we chose a unit slope). When mutation occurs, 
a mutant baby changes its x value to x±l with probabili- 
ties Af^(a;) = (l±^^(a;))/i/2, respectively. Random func- 
tions ^r and ^p are assumed to be uniformly distributed 
between ± A^.^ and they are "frozen" , or quenched, for 
a given realization. 

An exact description of this Markov process involves a 
very large (strictly speaking, infinite) set of master equa- 
tions for the probabilities of all possible population con- 
figurations. Significant simplification can be achieved in 
the limit of small mutation rate /i when all population is 
highly localized near (x) and only two neighboring sites 
x,x + 1 are usually occupied simultaneously (the proba- 
bility to occupy simultaneously 3 or more sites is smaller 
by a factor of /i) . The number of master equations then is 
reduced to N (see 0). It turns out that a much simpler 
system of only 2 individuals exhibits similar non-trivial 
features as A'^-individual model. This system in the limit 
of small /i can be described by only two equations for 
probability / to find a particle at site x, and g, the prob- 
ability of finding one particle at site x and another at 
x + 1, 

dtfix) - -^li?(a;)/(x) + ^R{x) {g{x) + g{x - 1)) , (1) 

dtg{x) = ^ (i?(a;)A/+(a;)/(x)-f R{x + 1) x 

xM-(a; + l)/(x + 1)) - ^ {K^) + R{^ + 1)) 9{^)- (2) 

The first term in r.h.s. of dl^) reflects the probability of 
mutation events leading to transition from the collapsed 
state when both particles occupy site x to either of dis- 
tributed states (x, X -\- V) ox (x — 1, x). The second term 
describes a reverse process of replication at site x which is 
followed by elimination of a particle at a; ± 1 and collapse 
of the population to the site x. Terms in the right-hand 
side of Eq.(p|) have the same origin. For ^ ^ 1, dtgix) is 



small in the asymptotic limit, and g is enslaved to /, 

_ R{x)M^{x)j{x) + R{x + \)M-{x + 1)/(.T + 1) 
^^^^ i?(.T) + i?(x+l) 

. Substituting g(x) into Eq.(Q), taking a continuous x 
limit, and keeping only linear in ^ terms, we obtain a 
single equation for /, 



dt!{x) 



3^9:^ [(1 



2a;C^ + d^ir)!{x)\ 



1 



(3) 



+ f5|[(x + e. + ^U/(^ 

In a similar fashion but with more cumbersome algebra, 
an equation describing the probability distribution /(x, t) 
for a population of N individuals to be localized at point 
X at time t can be derived in a small /i limit. In the limit 
of TV > 1 it reads 

dt!{x) = -|a, [(TV + 2.<^ + Nd^^r)f{x)] 

+ ^dl[{x + Cr + ^^,)f{x)]. (4) 

For Ar.p = this equation coincides with one derived 
in Ref. B and describes the increase of average fitness 
of N individuals as a result of random mutations and 
selective pressure. In the limit of smooth landscape the 
average fitness (x) = J xf{x)dx grows linearly in time 
with the rate Vb = ^ 0. A general solution to Eq.(||) 
in the presence of random fluctuations is not available, 
however one can estimate (x) using general methods of 
stochastic kinetics in random media [|lO| . To this end, 
note that Eq. (JSJ) has a form of Fokker-Plank equation for 
the dynamics of a single particle in the fitness space un- 
der the global bias TV/i/2 in the presence of the random 
quenched potential. Note that in our case the temper- 
ature itself depends on x. The corresponding dynamic 
process is governed by the Langevin equation: 

dx N 

— ^-(l + d^^r) + x^,. + vix,t), (5) 

where stochastic term r](x, t) has the following correlator 
{ri{x,t)ri{x,t')) = 2T(5(t — t') with local "evolution tem- 
perature" T = ^{x + £,r + N^f^/2), and we rescaled time 
lit -^ t. If the stochastic term were absent, the particle 
would have been pinned by local minima of the quenched 
potential. However, in the presence of "thermal fiuctu- 
ations" ri{x,t) the system will evolve via a sequence of 
thermally activated jumps from one minimum to a neigh- 
boring one favored by the global positive bias (selective 
pressure). This type of motion is well known in the dy- 
namics of disordered media and is usually referred to as 
creep g. 

We discuss the cases £,fj. ^ 0,£,r — and ^^ = 0, ^r ^ 
separately, since they exhibit different asymptotic behav- 
ior. Let us first focus on the case ^^ 7^ 0, ^^ = 0. It gener- 
alizes a well-known Sinai diffusion problem for the parti- 
cle subject to a random force field. The particle displace- 
ment obeys ultra-slow sub-diffusion law (x) r^ (logt)^ for 



zero average driving |ri|] , and the power-law dependence 
(x) ~ f^ in the driven case |12|. Our situation is more 
compUcated since the evolution temperature and the 
magnitude of disorder depend on x explicitly. To obtain 
the time-dependence of (x) we generalize the approach of 
Ref. Q . According to the general theory of stochastic 
growth fin] the evolution rate V — {x) oc f{U[x]), where 
U is the typical potential barrier controlling the evolu- 
tion, and / is the corresponding quasi-equilibrium distri- 
bution function, given by (we consider the limit x ':§> N 
and therefore drop the random term in the expression for 
evolution temperature) : 



f{x) ^ exp 2 / U^')dx' + {N-1) log(a:) 



(6) 



rM2\ 



Taking into account that the typical height of the po- 
tential barrier corresponding to the displacement x is 

'x/3 (by assumption £,m is 



uniformly distributed between ±A^) one finds the fol- 
lowing estimate for the typical time to evolve over the 
distance x. 



it) ^ V-' ^ l/fiUM) 

- exp{2Af,^/^-{N -l)\ogx). 



(7) 



Inverting Eq. (H) we obtain the following evolution law 
(x) ~ -^ (log(i + Co) + iN-l) log(x))' (8) 

where the constant Cq is determined by the initial con- 
dition for (x) at i = 0. In the limit {x) —^ oo, the final 
asymptotic behavior of Eq. m) coincides with the Sinai 
diffusion law ||ll|] (x) '^ ^(logi)^. For large N however 

this regime realizes at enormously large times t and dis- 
tances (x) {{x) ^ N'^ or t ^ e^), and can be hardly 
observed in numerical experiment. In the intermediate 
asymptotic regime (1 -C {x) ^ N'^), the last term in 
Eq. (ph is a leading one, and taking into account that 
log(a;) is a slow function as compared to (x), we obtain 
the following evolution law: 



N 

(x) ~ -rj log(^t + Co)+Xo, 



(9) 



(we restored original time units). Constants C'o and xq 
themselves depend on N and A^, and a crossover to 
disorder-free (A^ — > 0) linear growth of the fitness is re- 
covered by taking the limit fit <^ Cq and Co ^ oo. Our 
numerical results for discrete model simulations confirm 
analytical formula (B). Fig. 1 illustrates smooth loga- 
rithmic dependence {{x{t))) which was obtained numer- 
ically for N = 50, fj, = 0.02, Ayj = 0.7 averaged over 
50 statistically-independent configurations of quenched 
disorder (double angular brackets indicate computing a 



mean fitness of the entire population and ensemble av- 
eraging over different realizations of ^). As expected, 
for large t the rate of the evolution slows down which 
corresponds to the aging phenomenon known in glassy 
dynamics |0]. In the Fig. 1 we also show the dynamics 
of the mean fitness of the population (x) for a single run. 
This dependence is characterized by long periods of rel- 
atively steady fitness level (stasis) interrupted by abrupt 
jumps of the fitness. During these long intervals of stasis 
the population is trapped in regions of low probability 
of escape towards larger x (large negative ^^). Such an 
inhomogeneous pace of evolutionary changes is known in 
biology as a punctuated equilibrium hypothesis (see ^). 
We performed numerical experiments with our model for 
various N,fi, and A^. Fitting the results to the loga- 
rithmic dependence {(x)) = A\og{t + to) + xo yields a 
good agreement with formula (|9|). Figure 2 shows the 
dependences of AA'^ on N for several sets of parameters 
N, A^. All the graphs collapse close to a single straight 
line with slope 1 in log-log coordinates, as predicted by 
the theory. 

Now we turn to the case ^^^ = but ^r 7^ where par- 
ticle dynamics is governed by the diffusion in the random 
potential ^r rather than the random force field x^^. For 
a statistically uniform system (x-independent tempera- 
ture and magnitude of the quenched disorder) with the 
Gaussian statistics of the disorder ^r, the evolution rate 
V to the leading order is given by the following expres- 
sion (compare with the thermally activated hopping rate 
from 10) 



V ~ exp (-D/t^ 



(10) 



where D is the variance of the random potential and r is 
the temperature. In the present context both local vari- 
ance D and the "evolution temperature" r depend on x, 
however we expect that this formula holds adiabatically 
if the typical waiting time V^^ is large as compared to 
the relaxation time Vq" . The average evolution tempera- 
ture (r) = nx/2, and variance is D = V(^{0) = FffA^/S. 
The evolution rate is thus determined by the following 
asymptotic expression 



V '^ Vq exp 



W^ Al 
3fi^ x"^ 



7V2 A^' 

^o'^-pI-^^I (11) 



(here we added the pre-exponential factor Vq in order to 
match the results with the disorder- free case). Accord- 
ing to Eq.(O), the evolution rate V grows monotonously 
with X and finally approaches the disorder-free rate Vq- 
Therefore, in a long-time limit the quenched potential ^r 
becomes irrelevant because the growing with x evolution 
temperature smears out the pinning potential and further 
promote the growth of V . However, for x not too large 
{x <ti NAr) one can distinguish two different regimes 
of the (initial) evolution. If A^ 3> 1, in the absence 
of the thermal fluctuations a particle would be pinned 



by the random potential. Due to thermal fluctuations 
the particle will creep toward larger x, slowly increasing 
the velocity from y = 0. In contrast, if Aj. ^ 0(1), 
i.e. above the depinning threshold, the particle starts to 
move with the finite velocity V" < Vo, and after some 
time achieves asymptotic velocity Vq. Thus, the model 
exhibits "thermal depinning" and ignores completely fit- 
ness fluctuations for large average fitness values. Both 
regimes are seen in numerical simulations of the discrete 
model. Shown in Fig. 3 are {{x{t))) for different val- 
ues of Ar. We see that for large t all lines have the same 
slope which coincides with the evolution rate on a smooth 
landscape Vb = N^/2 (see also inset in Fig. 3). Inter- 
mediate regime exhibits a wide range of evolution rates 
which depend on A^. 

In conclusion, we have shown that two types of 
quenched randomness in the fitness space have a differ- 
ent effect on the population evolution. Quenched ran- 
domness of the replication rate slows down the evolution 
only at the initial stage since the increase of the "evo- 
lution temperature" eventually smooths out ruggedness 
and evolution proceeds at a rate corresponding to the 
smooth landscape case. However, in the long-time limit 
quenched disorder of the mutation rate being amplified 
by the increasing replication rate, dominate the evolution 
temperature growth. In the limit x, i — > oo the evolu- 
tion exhibits ultra-slow logarithmic growth of averaged 
mean fitness {{x)) . Individual runs are characterized 
by long intervals of almost constant mean fitness (x) in- 
terrupted by spontaneous changes, a dynamics which is 
usually interpreted as punctuated equilibrium. In this 
work, we have made simplest ad hoc assumptions re- 
garding the statistical properties of the fitness landscape. 
An important task for a future work will be to estab- 
lish connections between "standard" multi-dimensional 
sequence space evolution models such as NK-model and 
one-dimensional fitness space model described here. 
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FIG. 1. Averaged over 50 runs mean fitness {{x{t))) for 
iV = 50 and ^ = 0.02, A^ = 0.7, and A^ = (solid line). 
Dashed line shows {x{t)) for a single run. 




FIG. 2. Scaling of pre-Iogarithm factor A with N and A;j, 
obtained by fitting {{x{t))) — Alog{t + Co) + xq. Solid line 
indicates a unit slope. 
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FIG. 3. Averaged over 100 runs mean fitness ((x(t))) for 
A*' = 100, /i = 0.1, Ap = and different Ar. Inset shows 

V = a((^». 



